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1. 


MOTIVATION 


A  family  of  simulation  models  has  been  developed  to 
provide  Monte  Carlo  simulation  of  atmospheric  variables 
in  time  and  space.  These  models  are  based  on  the 
superimposition  of  waves  to  approximate  a  standardized 
Gaussian  field.  Three  parts  of  these  models  must  function 
accurately  if  the  generated  fields  are  to  imitate  the 
atmosphere : 

a.  Convergence  to  a  normal  distribution  of 
the  superimposed  waves. 

b.  Transformation  of  the  normal  distribution 
to  the  variables'  climatological  distribu¬ 
tion. 

c.  Correlation  structure  in  the  generated 
Gaussian  field  must  closely  approximate 
the  variables'  equivalent  normal  correla¬ 
tion  structure  in  time  and  space. 

It  is  the  third  requirement  that  this  report  addresses. 
Furthermore,  the  results  are  applicable  to  other 
statistical  atmospheric  problems  including  objective 
analysis,  inferential  testing,  coverage  distributions, 
and  extremes. 
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2. 


WAVE  GENERATORS 


Here,  a  wave  is  defined  as  a  periodic  function  in 
one  or  more  dimensions: 

Wave  (u)  =  Wave  (u  +  L)  (1) 

where  wave  is  some  arbitrary  waveform,  e.g.,  a  cosine  or 
a  sawtooth,  L  is  the  wavelength  and  (1)  is  true  for  any 
value  of  u.  Initial  consideration  is  limited  to  plane 
waves,  that  is,  u  is  a  linear  function. 

In  3  dimensions;  u  =  CxX  +  Cyy  +  CzZ  +  Uo  (2) 

In  D  dimensions:  u  =  c*x  +  Uo  (3) 

where  the  c's  are  direction  cosines  of  a  line  (the  cosine 
of  the  angle  between  each  axis  and  the  line)  with 
positive/negative  orientation  through  the  origin,  c  x  is 
a  dot  product  between  c,  a  row  vector  of  direction 
cosines,  and  x,  a  column  vector,  and  Uo,  0  <  Uo  <  L,  gives 
the  phase  at  the  origin.  Equation  (3)  defines  a  plane 
perpendicular  to  the  line  specified  by  the  direction 
cosines.  As  illustrated  in  Figure  1,  the  value  of  the 
wave  is  the  same  everywhere  on  this  plane. 

Equation  (1)  along  with  (2)  or  (3)  define  the  value, 
V,  of  the  wave  everywhere  in  the  space  considered.  If  a 
set  of  such  waves,  each  with  different  direction  cosines 
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Figure  1.  Geometry  for  a  three-dimensional  plane  wave. 
Ray  along  line  from  0  to  P  is  defined  by  the  direction 
cosines  of  the  angles  Cj,  C2,  and  C3.  The  wave  has  the 
same  value  everywhere  on  the  plane  perpendicular  to  OP 
and  passing  through  A,  B,  and  C. 


and  origin  phases,  is  added  together,  an  elaborate 
pattern  results.  If  the  direction  cosines  and  origin 
phases  are  randomly  chosen,  the  pattern  will  be  random, 
but  will  have  an  underlying  correlation  structure 
depending  on  the  type  of  wave  and  method  of  selecting 
direction  cosines  and  origin  phases.  Also,  as  the  niomber 
of  waves  goes  to  infinity,  the  resultant  field  may  have 
fractal  properties. 
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3.  CORRELATION  STRUCTURE  FORMULAS 

To  simplify  the  formulas,  waves  are  first 
standardized  so  that  over  a  wavelength  the  mean  is  zero 
and  standard  deviation  is  one.  When  m  (m  ^  2)  random- 
independent  wave  fields  are  added  together  with  weights 

{Wi)^2^  the  resultant  field,  F, 

F^c)  =  Z  Wavej  (cj  ’  ^  (4) 

i  =  1 

will  have  a  variance  of  one  if. 


m 

Z  W.  =1 
i  =  l 


(5) 


and  the  waves  have  been  standardized.  Variance  can  be 
defined  in  two  (at  least)  ways  for  the  field  F.  One  way 
is  to  use  all  values  of  x  in  a  single  realization  to 
calculate  the  variance.  The  other  way  is  to  calculate 
variance  at  a  single  point  using  multiple  realizations. 
Equation  (5)  is  true  for  both  cases. 


Given  multiple  realizations,  suppose  the  correlation 
between  the  values  of  wave  i  of  m  waves  at  Xi  and  Xa  is  ri 
(xi,X2)/  and  that  the  i  waves  are  independent  of  each 
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other,  then  the  correlation,  rp,  between  values  of  F  at  Xi 
and  X2  is: 


m 


rp=I  Wj  r 
i  =  l 


(6) 


Proof:  The  correlation  between  a  set  of  standardized 

variables  v(xi)  and  v(x2)  is: 


K 

Z  v(x,)  v(x,) 
k  =  l 


(7) 


where  the  summation  is  over  multiple,  1  to  K,  realiza¬ 
tions.  Combining  (4)  and  (7)  : 

Km  m 

rF(x,,x,)=Z  [Z  (x,)][Z  Wl'^v  (X3)] 

k  =  l  i=l  ‘  ‘  j  =  l  J  J 

All  cross  terms  (VjVj,  I  j )  will  go  to  zero  as  K  goes  to 
infinity  since  the  waves  are  independent  of  each  other. 
Thus,  with  only  J=i  terms  (8)  becomes: 


K 

m 

h  (x>,  X,) 

=  Z 
k  =  l 

Z 

i  =  1 

(x,)][WP  Vj(x.)] 

(9) 

m 

K 

=  ZW; 
i  =  l 

Z 

k  =  l 

Vj  (X|)  V.  (X,) 

(10) 

and  using  (7) 

becomes 

(6)  . 

The  importance 

of 

this  result  is  that 

many 

waveforms,  each  with 

its 

own  correlation  structure 

,  can 
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The  importance  of  this  result  is  that  many 
waveforms,  each  with  its  own  correlation  structure,  can 
be  superimposed  to  provide  a  great  variety  of  correlation 
structures.  Further,  since  the  correlation  structure  is 
generated  from  a  set  of  realizations  (ignoring  the 
trivial  case  v-0  everywhere) ,  the  correlation  structure 
resulting  from  a  superimposition  of  waves  will  have  a 
positive  definite  matrix  of  correlations  for  any  set  of 
points.  Thus,  any  of  the  resulting  correlation  structures 
are  eligible  for  consideration  in  objective  analysis. 
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4. 


HOMOGENEOUS  ISOTROPIC  CORRELATION 


If  the  correlation  depends  only  on  the  distance 
between  two  points,  not  on  where  or  in  what  direction 
distance  is  measured,  the  correlation  structure  is  said 
to  be  homogeneous  and  isotropic.  Such  a  structure  results 
when  a  field  is  the  result  of  many  waves  superimposed 
with  their  directions  uniformly  random,  that  is,  any 
direction  is  as  equally  likely  as  any  other. 

One  way  to  pick  a  random  direction  is  to  randomly 
pick  a  point  in  a  sphere  and  to  use  the  ray  from  the 
center  of  the  sphere  to  the  point  to  specify  a  direction. 
■A  random  point  in  a  sphere  can  be  chosen  by  picking  a 
random  point  in  a  cube  enclosing  the  sphere  using  uniform 
random  variables,  but  only  using  the  point  if  it  lies 
within  the  sphere.  This  algorithm  is  described  in  detail 
in  Gringorten  and  Boehm  (1987).  This  and  other  algorithms 
are  discussed  in  Marsaglia  (1972) . 

4.1.  Direction  Cosines 

The  direction  cosines  specify  the  direction  of  a 
linear  wave.  Given  that  any  direction  is  equally 
probable,  what  is  the  probability  density  of  the 
direction  cosines? 
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Given  the  area  of  a  sphere  (in  2  dimensions  "area" 
is  the  length  of  arc  on  a  circle  while  in  D  dimensions 
"area"  is  the  boundary  of  a  hypersphere)  less  than  a 
given  direction  cosine#  the  change  in  area  per  change  in 
direction  cosine  gives  the  probability  density#  fco#  that 
a  specific  direction  cosine#  c#  would  be  chosen: 


1 

step  function  .5  at  -1 

spike  at  -1  and  1 

(11) 

2 

[sin”^(c)/n]  +  1/2 

1/  [n(l  - 

(12) 

3 

(c  +  l)/2 

1/2 

(13) 

4 

sin"Mc)+c(l-c^)^^^  +  1/2 
n 

2(1-  c^) 
n 

(14) 

5 

2  +  3c  -  c^ 

3(1  -  c^)/4 

(15) 

n 


^c+l  (  ^  (1  -  c^) 

—  -  (16) 

B(1/2,(D  - 1)/2 

where  Ix(a#b)  is  the  incomplete  beta  function  using  the 
notation  of  Eq  (26.5.1)  of  Abramowitz  and  Stegon  (1964)# 
and  B  is  the  beta  function: 

B[l/2#  (D-l)/2]  =  r(D/2) 

r(l/2)r( (D-1) /2) 

The  left  side  of  (11)  to  (16)  is  the  integral  of  the 
right  side.  In  the  general  case#  (16)#  the  form  of  the 
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incomplete  beta  function  is  obtained  from  the  right  side 
by  noting  that, 

(1_C)  [(d-3)/2)  (i^C)  ^  ((d-3)/2] 

and  changing  the  variable  of  integration  c=2x-l. 

The  direction  cosine  density  is  symmetric  and 
therefore  the  odd  moments  are  zero  and  the  raw  moments 
are  equal  to  moments  about  the  mean.  Given  the  form  of  a 
beta  distribution  through  the  change  of  variable  s=l-c^, 
moments  of  the  direction  cosine  density  are  easy  because 
the  moment  integral  itself  is  given  by  a  beta  function. 
The  mth,  moment,  p'"’’,  is: 

^  B[(m  +  l)/2,(D-l)/2]  (^9) 

^  B[l/2,(D-l)/2] 

As  D  becomes  large,  the  distribution  of  the 
direction  cosine  approaches  a  normal  distribution  with  a 
mean  of  zero  and  variance  of  1/D. 

To  obtain  values  along  a  single  axis  in  D  space,  it 
is  not  necessary  to  generate  direction  cosines  for  all 
the  axes.  Since  the  correlation  is  isotropic,  any  line 
can  arbitrarily  be  aligned  as  the  first  axis  then  all 
terms  in  (2)  or  (3)  except  the  phase  at  the  origin,  Uo, 
and  the  first  direction  cosine  term  are  zero.  In  this 
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case,  the  cimulative  probabilities,  the  left  side  of  (11) 
through  (16),  can  be  used  to  generate  a  random  direction 
cosine  using  a  uniform  random  number: 

To  directly  generate  values  along  a  second  or  third 


axis. 

the 

conditional 

probability  must 

be 

used. 

The 

conditional 

density  is  the 

probability 

density  of 

the 

next 

lower 

dimension 

times 

a  factor  that 

reduces 

its 

range. 


=  (1  (20) 

D 

=  -  z  (21) 

j  =  m  +  l  •’ 

4.2.  Isotropic  Correlation  in  D  Dimensions 

Consider  the  correlation  between  values  at  two 
points  separated  by  a  distance,  <5,  along  a  one 
dimensional  wave: 

rj(^)  =  {  j*  Wave(x)  Wave(x+^)  d\-ju^}/a 
0 

where  x  and  5  are  in  units  of  wavelength,  fj  is  the  mean 
of  the  wave,  and  6  the  standard  deviation.  Two  ways  of 
viewing  (22)  are  equally  valid.  First,  consider  the  two 
points  as  fixed  and  x  as  the  phase,  i.e.,  integrating 
over  all  phases  zero  to  one.  The  second  way  is  to 
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consider  the  phase  fixed  and  x  as  a  location,  that  is, 
integrating  overall  locations  from  zero  out  to  one 
wavelength. 

If  multiple  waves  are  used  as  in  (4),  then  the 
resulting  correlation  functions  can  be  calculated  by  (6) . 
If  the  probability  density  of  a  given  wave  is  used  for 
the  weight,  the  resulting  correlation  function  is  an 
expected  value,  e.g.,  the  asymptotic  result  of  many 
realizations . 

Further,  if  the  probability  density  used  is  the 
probability  density  of  a  direction  cosine  as  in  (11)  to 
(16)  above,  the  resulting  generality  separation  distance, 
5,  may  be  measured  along  the  first  axis  and  the  first 
direction  cosine  completely  defines  the  equivalent 
wavelength,  1/c,  along  that  axis. 

Correlation  as  a  function  of  separation  distance,  5, 
can  be  calculated  by  integrating  [vice  summing  as  in  (6)] 
over  all  possible  values  of  the  direction  cosine  weighted 
by  that  value's  probability  density; 

1 

=  J  fcD(c)  dc  (23) 

-1 
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The  integral  is  sometimes  easier  to  integrate  if  the 
lower  limit  is  set  to  zero  and  the  result  multiplied  by  2 
(because  of  symmetry)  and  the  distance,  s,  s=<5/c,  is 
substituted  as  a  change  of  variable: 

.  1 

^  T  I  ri(s)  ds  (24) 

*  0 

(24)  along  with  (11)  through  (16)  and  (22)  were  used 
to  find  correlation  functions  for  several  waveforms. 


4.3.  Correlation  of  Sine  Waves 

Sine  and  cosine  waves  have  a  mean  of  zero  and 
variance  of  1/2.  The  wavelength  of  sine  (and  cosine) 
waves  are  generally  expressed  in  radians  so  that  if  <5  is 
in  wavelengths,  arguments  must  be  multiplied  by  2n.  The 
one-dimensional  correlation  is  simply: 

risin(^)  =  cos(2;r<5) 

from  which  the  higher  dimensional  correlations  can  be 
calculated: 


sin(2;r<J) 


2k5 


(26) 

(27) 
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Separation  Distance  a 


Sine/Cosine  Wave  Correlation 

Figure  2.  Sine  or  cosine  wave  correlation  as  a 
function  of  separations  distance  for  waves  with  uniform 
distributions  of  direction  and  phase  in  2,  3,  4,  and  5 


dimensions . 


r4sin(<^)  =  2  l,(27,5)l{ln5) 


(28) 


3sin(2;j5)  -  2k5co^z5)  (29) 

TssinW  = 

where  Jn  is  the  nth  order  Bessel  function  of  the  first 
kind  (Abramowitz  and  Stegon,  1964) .  Bessel  functions  can 
not  be  evaluated  directly,  however,  series  expansions, 
approximations,  and  computer  algorithms  are  available. 

Since  other  waveforms  can  usually  be  represented  as 
a  series  of  sine  and  cosine  waves  via  the  Fourier 
transform,  the  correlation  functions  can  also  be 
represented  as  a  series  of  sine  correlation  functions. 
However,  integrating  the  forms  directly  most  often 
results  in  a  simpler  form. 


4.4.  Correlation  of  Savtooth  Waves 

The  sawtooth  wave  is  discontinuous  at  integer  values 
of  wavelength,  <5,  so  that  it  is  convenient  to  define  f  as 
fraction  of  wavelength  (or  phase) : 


f  =  ^  -  int(^) 

where  int{(5)  is  the  largest  integer  in  5.  Note:  int  is 
defined  this  way  in  Basic  but  only  for  positive  values  in 
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FORTRAN.  The  value  of  the  sawtooth  wave  is  equal  to  the 
phase, 

Saw(<^)  =  f  (31) 

Integrals  are  generally  evaluated  piecewise  between 
integer  values,  so  that  when  the  integral, 


JSaw(x)Saw(x  +  <5)  dx  =  J  x(x  +  ^)  dx  +  J  x(x+^-l)  dx 
0  0  1-^ 

=  (2  -  3^  +  35^)16 


(32) 


is  combined  with  the  mean  of  1/2  and  variance  of  1/12  for 
the  sawtooth  wave  the  correlation  is  given  as  a  function 
of  separation  distance: 


r,saw(^)  =  1-6^  +  6S\ 

or, 

x,^i5)  =  1  -  6f  +  6f^ 
Integration  of  (24) 
and  gives  in  general: 


S  <  1 


(33) 


f  =  ^  -  int(<J) 

using  (33)  is  straightforward 


*Dsaw  i  ^  —  1)  —  1 


_ ^  ^  6  o2 

(D-l)B(l/2,D/2-l/2)  D 


(35) 


In  fitting  models  to  observed  correlations,  the 
inverse  equation,  i.e.,  distance  given  correlation,  is 
sometimes  needed.  This  distance  is  unambiguous  only  as 
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far  as  the  first  minimum.  The  lower  dimensional 
equations,  their  inverses  and  first  minimums  (rll^*"  min) 
are : 


r2saw(^l<^<l)  =  l  -  —S  +3S^ 

n 

(36) 

5^i.x\5<2ln-)  =  Hk  -  [(r-l)/3  +  4/;r^f 

(37) 

r33.„(^|^<l)  =  1-3^ 

(38) 

^3saw(r|^<3/4)  =  3/4  -  (r/2  +  1/16)'^ 

(39) 

rteaw(^l^<l)  =  1-  —S  + 

n  1 

(40) 

^  -  [2(r  -  \)n  + 

37t  371  9n 

(41) 

For  5>  1,  values  of  correlation  can  be  obtained  by 
integrating  (24)  by  separating  6  into  two  parts, 

J  =  I  +  f,  where  I  =  int(^),  andf  =  S  -  int(<5),  (42) 

and  integrating  each  term  separately.  After  combining 
terms,  the  correlation  is  given  by: 

r^=  1-^  +  3^  +  6a  +  f)- 

^  (43) 

—  i  [i  sin-'(h  + 

TT  i=l  O 

=  (f  -3f'  +  2f)/S  (44) 
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4saw 


l.ls  +  -S^  +  6(I  +  f)  - 

TT  2 

-  i  [(2  +  ™’'4)i 

i=I  s  s 


(45) 


=  1  -  1.255  +  1.2  +6(1"  +I)-4.5IJ  - 

31(1"  +1.51  + 1/2)/<^  +  1(0.31^  +  0.751" +0.51"- 0.05)/ 

4.5.  Correlation  of  Triangular  Waves 

The  value  of  the  triangular  wave  can  be  written: 

Tri(J)=l-2f,  0<f<l/2 
=  2f  -  1,  l/2<f<l 

with  f  =  5  -  int  (5)  as  above.  The  convolution  integral  can 
then  be  written: 


1^  l/2-f 

J  Tri(x)Tri(x  +  f)  dx  =  J  (l-2x)[l-2(x+f)]  dx 
0  0 

1/2  1-f 

+  J  (l-2x)[2(x  +  f)-l]  dx  +  J(2x-l)[2(x+f)-l]  dx 

l/2.f  1/2 

1 

+  J  (2x-l)[l-2(x-l-f)]dx,  forf<l/2 

1-f 


=  1/3-2^"+-^",  for  f<  1/2 

3 


(48) 


(49) 
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Separation  Distance  a 


Sawtooth  Wave  Correlation 


Figure  3.  Sawtooth  wave  correlation  as  a 
function  of  separations  distance  for  waves  with  uniform 
distributions  of  direction  and  phase  in  2,  3,  4,  and  5 


dimensions . 


A  similar  set  of  integrals  for  1/2  <  f  <  results  in: 


jTri(x)Tri(x  +  <J)  dx  =  l-4f  + 

forl/2<f<l 

When  these  results  are  combined  with  the  triangular 
wave's  mean  of  1/2  and  variance  of  1/12,  the  correlation 
becomes : 


r,^(^|f<l/2)  =  1  -  24  f  +  32f^ 


r,^(<^|l/2<f  <1)  =  9  -  48f  +  72f^  -  32f' 


When  (50)  is  put  into  the  D-dimensional  correlation 
formula  (24),  the  following  are  obtained: 


r^i5\S<\l2)  =  1  -  US^  + 

371 


r3tri(^l^<l/2)  =  1  -  +  SS^ 


=  \-6s^  + 

IStt 
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For  separation  distances/  5/  greater  than  one-half/ 
the  integral  is  possible  but  becomes  increasingly 
complicated: 


r2,ri(<^)  =  S{sin-'(^)[(288i  +  72)/52  +  192i^+144p-6) 
i=0  O 

-  sin''  (^-^^)[(576i  +  288)^'  +  384  +  S16\^  +  288i  +  48] 

s 

+  sin''  (— )[(288i  +  216)^^  + 192  +  432  +  288i  +  54] 

s 

+  {5^ (128<J"  +  352  +  216i) 

+  [^^-(i  +  l/2)^]'^[-128^^-352(i2  +  i)-88]  (5g) 

+  -(i  +  l)^]''^(12&^^+352i^+488i  +  136)}  +R(f) 

where  R(f)  is  [Recall  from  (8)  f  =  5  -  int(5)  and  I  = 

int  (5)  ]  : 


R(f  I  f  >1/2)  =  sin  •'  (-)[(288I  +  72)  + 1 92  f  + 144  f  -  6] 

6 

.  sin-i  +  288)  <5^  +  384  f  +  576  f  +  2881  +  48] 

s 

+  )"'(128^"  +  352l'  +  216I)  +  2[S^ -(l  +  ^f][-ms^ 

-  352  (f  + 1)  -  88]  +  ;;r[(144I  + 108)^^  +  96  f  |  +  216 1"  + 1441  +  27] 


(57) 


R(f  I  f  <  1/2)  =  sin''  (4)[(288I  +  72)  <5^  + 1921'  + 144  f  -  6] 
o 

+  -l)^'^(\2ZS^  +  352f  +  2161)-  ;r[(144I+  36)  +  96f  +  72f  - 3] 


(58) 
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The  three-dimensional  triangular  wave  correlation 


function  is  easier  because  of  the  simpler  form  of 
three-dimensional  directional  cosine  density: 

r3^(^|f<l/2)  =(f-  8f^  +  SfyS 

r3^(^|f  >1/2)  =  -[(l-f)-8(l-f)^+  8(l-f)'']/^ 


the 

(59) 

(60) 
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1 


Separation  Distance  a 
Triangular  Wave  Correlation 


Figure  4.  Triangular  wave  correlation  as  a  function 
of  separations  distance  for  waves  with  uniform 
distributions  of  direction  and  phase  in  2,  3,  and  4 
dimensions . 
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5. 


DISCUSSION 


This  report  presents  the  method  of  calculating 
correlation  functions  for  periodic  functions  -  waves  -  in 
D-dimensional  space.  Explicit  formulas  for  correlation 
functions  of  several  important  waveforms  are  given.  These 
waveforms  include  the  sine/cosine,  the  sawtooth  wave,  and 
the  triangular  wave. 

The  correlation  function  for  triangular  waves  is 
nearly  the  same  for  as  that  for  sine  waves.  When  the 
graph  for  one  is  placed  on  top  of  graph  for  the  other, 
the  graphs  appear  to  be  identical.  This  similarity  has 
two  important  applications.  In  simulating  a  field  that 
has  a  sine  type  of  correlation,  the  triangular  wave  can 
be  used;  the  triangular  wave  being  an  order  of  magnitude 
faster  to  simulate  than  a  sine  wave. 

The  second  application  is  in  approximating  the  first 
order  Bessel  function.  The  2D  sawtooth  correlation  is 
particularly  applicable  in  inverse  Bessel  function. 

Often  objective  analysis  requires  correlation 
functions  of  a  certain  shape.  The  sawtooth  has  a 

correlation  similar  to  exp (-6^)  .  By  mixing  the  sawtooth 
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and  triangular  wave  correlation  curves  and  by  also  mixing 
the  wavelengths,  a  variety  of  correlations  is  possible. 
These  correlation  curves  are  guaranteed  to  regression 
using  correlations  will  always  (in  theory)  be  able  to  be 
inverted.  Further,  objective  analysis  is  mainly  concerned 

with  close  (less  than  one  5)  points.  For  these  points, 
the  resulting  mixture  is  a  simple  quadratic  polynomial 
that  is  very  fast  to  calculate. 

At  first,  sawtooth  correlation  functions  were  only 
available  as  the  result  of  numerous  simulations.  When 
these  were  to  be  combined  to  approximate  another 
correlation  function,  niomerous  simulation  runs  were 
required.  In  order  to  see  the  effect  of  varying  different 
weights  and  wavelengths,  graphs  from  these  different  runs 
were  placed  side  to  side  and  top  to  bottom.  They  filled 
the  whole  floor  of  a  room.  What  a  help  the  equations  in 
this  report  could  have  been  then. 


6. 
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